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We show that, when a single relaxation time lattice Boltzmann algorithm is used to solve the 
hydrodynamic equations of a binary fluid for which the two components have different viscosities, 
strong spurious velocities in the steady state lead to incorrect results for the equilibrium contact 
angle. We identify the origins of these spurious currents, and demonstrate how the results can 
be greatly improved by using a lattice Boltzmann method based on a multiple-relaxation-time 
algorithm. By considering capillary filling we describe the dependence of the advancing contact 
angle on the interface velocity. 
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I. INTRODUCTION 

There are many different physical systems for which 
contact line dynamics plays an important role. These 
range from large scale industrial processes, such as oil re- 
covery and coating, to micron scale phenomena, such as 
flow in microfluidic channels or the movement of drops 
across micropattcrncd surfaces. Despite its importance, 
contact line dynamics is still a challenging problem that 
has not been fully solved. This is partly because it is an 
inherently multiscale problem. The hydrodynamic flow 
of the fluids far away from the three phase contact line 
has to be matched to the microscopic details of the con- 
tact line motion In particular, a mechanism is needed 
by which the usual hydrodynamic, no-slip boundary con- 
dition can be violated in the vicinity of the contact line 

Si. 

One solution to this problem has been to consider a 
class of hydrodynamic models, called diffuse interface 
models [1, H, [1] , where the interface at the contact line 
has a finite width. As a diffuse interface is pushed across 
a surface it is pulled out of equilibrium. This leads to 
diffusive transport of fluid across the interface and hence 
an effective slip at the contact line. Here we shall con- 
sider a diffuse interface model for a binary fluid, where 
the equilibrium properties of the fluid are described by a 
Landau free energy functional. The equilibrium contact 
angle is controlled by a surface term in the free energy. 

The hydrodynamic equations of motion of diffuse in- 
terface models can be solved in many different ways, 
but in this paper, we shall focus on one particular 
mesoscopic modelling technique, the lattice Boltzmann 
method 0, H, 0- Lattice Boltzmann simulations have 
been succesfully used to study several contact line prob- 
lems. Examples include the spreading of drops on chem- 
ically patterned [l^l and superhydrophobic surfaces [ll|, 
modelling of contact angle hysteresis 12 1, and capillary 
filling in microfluidic channels [H, [ij, 15 1 . 

The majority of previous work on binary fluids has con- 
centrated on the case where the two components have 
equal viscosity. In this case a simple and widely used 
lattice Boltzmann approach, the BGK algorithm, works 



well, agreeing with analytical results for the equilibrium 
contact angle. However, if the two binary components 
have different viscosities, this is no longer the case: we 
shall show that the equilibrium contact angle is predicted 
incorrectly by the algorithm because there are severe spu- 
rious velocities in the steady state. 

There are several important aspects of the behaviour 
of multicomponent fluids where it is essential or highly 
desirable to be able to model a two component fluid, 
where the components have different viscosities. These 
include drops moving on surfaces, where the viscosity of 
the surrounding fluid must be substantially smaller than 
that of the drop to access a rolling regime, instabilities 
such as viscous fingering which are driven by a viscosity 
difference between the two fluid components, and cap- 
illary flUing, where the simple theories assume that the 
displaced fluid has zero viscosity. 

In this paper we identify two primary reasons for the 
spurious currents in BGK lattice Boltzmann simulations 
of contact line hydrodynamics. We show how the spuri- 
ous effects can be greatly suppressed by using an algo- 
rithm based on a multiplc-rclaxation-time lattice Boltz- 
mann approach. We demonstrate that the new algorithm 
gives excellent agreement with theory, both for the equi- 
librium contact angle, and for the advancing contact an- 
gle, measured in capillary filling simulations. 

The paper is organised in the following way: We be- 
gin, in Sec. |TT1 by introducing the free energy model for a 
binary fluid system. We summarise two different lattice 
Boltzmann implementations that can be used to solve 
the binary fluid's hydrodynamic equations of motion, the 
BGK model and a multiple-relaxation-time method, in 
Sees. IIIII A and B. In Sec. IIVI we measure the equilib- 
rium contact angle for a drop on a surface and flnd that, 
for the BGK approach, it deviates from the value pre- 
dicted theoretically if there is a viscosity ratio between 
the two phases. This deviation is caused by anomalous 
spurious currents near to contact points. In Sec. |Vl we 
discuss the origin of the spurious currents: flrstly long 
ranged effects and secondly non-zero velocities induced 
by the bounce-back boundary conditions. We propose, in 
Scc. lVIi an algorithm based on a multiplc-rclaxation-time 
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lattice Boltzmann implementation [T^, [l^ which signif- 
icantly suppresses the spurious currents at the contact 
point. In Sec. IVIll we use simulations of capillary filling 
in smooth channels to measure the dependence of the 
advancing contact angle on capillary number. Finally, in 
Sec. IVIIH we summarise the results and conclude. 



III. LATTICE BOLTZMANN 
IMPLEMENTATIONS 

The equations of motion ([6|)-(l8]) can be solved using 
a lattice Boltzmann algorithm. We shall consider two 
different lattice Boltzmann approaches in this paper and, 
for convenience, we summarise both here. 



II. MODELLING THE BINARY FLUID 

The equilibrium properties of a binary fluid can be 
described by a Landau free energy functional 



F 



{i' + ^\v4>f}dv- 
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where the bulk free energy density ip is taken to have the 
form 



?A = ^p\np + a (- 



1 wi,2 



(2) 



The first integral in Eq. ([T]) is taken over the volume of 
the system, p is the fluid density and 4> is the concentra- 
tion, c is a lattice velocity parameter, described below 
and a is a constant. This choice of t/j gives binary phase 
separation into two phases with = ±1. The k term in 
Eq. ID) represents an energy contribution from gradients 
in (j) and is related to the surface tension between the two 
phases through 7 — ^8Ka/9 0]. 

The second integral in Eq. ([1]) is over the system's sur- 
face and is used to model the fluid-solid surface tensions 
p^ . The parameter ^ determines the contact angle. 

Taking the functional derivative of Eq. ((T]) with respect 
to gives the chemical potential 



(3) 



which is constant in equilibrium. Minimisation of the free 
energy also shows that the gradient in cf) at the boundary 
is = ^/k 0], where the partial derivative d± is 

taken in a direction normal to the surface. 

The dynamics describing how the fluid approaches 
equilibrium are determined by the pressure tensor 

= {po - - f |V0|2) 6o,p + Kdc,(t>dp(t), (4) 

where the bulk pressure is 

The hydrodynamic equations for the binary fluid are 



dtP + da{pVa) = 0, (6) 
dtipVp) + daipVaVfs) = -daPap (7) 

+da {vp {dpVa + daV/s)} , 



(8) 



where v is the fluid velocity, v is the kinematic viscos- 
ity and M is a mobility. The equilibrium properties of 
the fluid appear in the equations of motion through the 
pressure tensor and the chemical potential. 



A. Single-relaxation-time lattice Boltzmann 

To implement a lattice Boltzmann algorithm for a bi- 
nary fluid in two dimensions the system is divided up 
into a square grid of points, with two particle distribu- 
tion functions fi{r,t) and gi{r,t) on each point. The 
label i denotes a particular lattice velocity vector e^, 
defined by Gq = (0,0), ei,2 = (±c, 0), 63^4 = (0, ±c), 
65^6 = (±c, ±c), and 67^8 = (=Fc, ±c). The lattice veloc- 
ity parameter c is defined by c = Ax/ At, where Ax is the 
spacing between nearest neighbouring points and At is 
the time step. The physical variables are obtained from 
the particle distribution functions using 



(9) 



The time evolution equation for the particle distribu- 
tion functions, using the standard BGK approximation 
Q, can be broken down into two steps. The first is a 
collision step 

(10) 
(11) 



/;(r,t) = /,(r,i)-^[/,-/,n 
5Kr,t) = 5.(r,t)-^[5.-5r]- 



This is followed by a streaming step, which moves parti- 
cles along their corresponding lattice velocity vector di- 
rection 

/,(r + e,At,t + At) = /;(r,t), (12) 
g,{r + e,At,t + At) = g^(r,t). (13) 

/f^ and (7j^^ are equilibrium distribution functions, de- 
fined as a power series in the velocity. A summary of 
our choice of equilibria, which is motivated to help re- 
duce spurious velocities around interfaces p7| . is given 
in Appendix 1X1 Note, the inter-coupling between fi and 
gi comes through f^'' and g^'^. In particular, the large 
variation in </) at an interface influences fi by dacj) ^md 
terms in f^"^. 

A Chapman Enskog expansion can be performed to 
show that the lattice Boltzmann collision (fTUl [TT|) and 
streaming [T5)l operations lead to the hydrodynamic 
equations ([5])-® in the limit of long length and time 
scales 0. The relaxation parameters Tp and arc re- 
lated to the kinematic viscosity and mobility through 



= At^- (rp - i) 
I) 



M ^ Atr (r0 



(14) 
(15) 



where F is a tunable parameter that appears in the equi- 
librium distribution. 
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B. Multiple-relaxation-time lattice Boltzmann 

We next summarise the basic methodology behind the 
multiple-relaxation-time (MKT) lattice Boltzmann ap- 
proach. More details are given in [l^, [l^. The idea be- 
hind multiple-relaxation-time lattice Boltzmann is that 
different relaxation parameters are used for different lin- 
ear combinations of the distribution functions. In partic- 
ular, the relaxation parameters responsible for generating 
the viscous terms in the Navier-Stokes equation ([7]) are 
set to Tp, those connected to conserved quantities to oo, 
and all others to 1. 

In multiple-relaxation-timc lattice Boltzmann the col- 
lision term [fi — f^''] on the right hand side of the 
lattice Boltzmann equation for fi (jTU]) is replaced by 



(16) 



where the particle distributions fi and f- are written as 
column vectors and M is a matrix. One possible choice 
for M is 
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This matrix performs a change of basis. The new basis 
is designed to contain more physically relevant variables. 
For instance, when the first row is dotted with f the 
density p ~ fi is obtained. Similarly, the fourth and 
sixth lines calculate the momentum densities pu^ and 
puy, respectively. Each of the rows in M is mutually 
orthogonal so the inverse follows easily as 



m; 



1 



M ■ 



(18) 



The matrix S in Eq. is diagonal and has the 

elements 



S = diag(0,l,l,0, 1,0, l,w,w) 



(19) 



where uj = l/^p now determines the fluid viscosity in 
Eq. Note that some of the values are zero. This 

choice is arbitrary as these modes correspond to the con- 
served moments; J^i ^Iji [fi ~ /f] = for j ~ 0, 3, 5. 
The choice of unity for the other terms will be justified 
later. 

It is not necessary to adopt a multi-relaxation ap- 
proach for gi as there is an independent parameter F, 
which can be varied to change the mobility of particles 
in Eq. (|15p . Therefore, even when using a MRT ap- 
proach, we set distribution gi to evolve according to Eq. 
(HH) with T0 = 1. 
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FIG. 1: (Color online) The simulation geometry. The contact 
angle is defined by 6"'' . The crosses denote the points which 
are fit to the arc of a circle to obtain a numerical estimate of 
the equilibrium contact angle 0°'. The grid is not to scale. 



IV. MEASURING THE EQUILIBRIUM 
CONTACT ANGLE 



In this section we check the extent to which the single 
relaxation time lattice Boltzmann method gives the cor- 
rect equilibrium contact angle at the contact point. We 
find that, for a fiuid with a constant viscosity modelled 
using Tp = Ta = Tp = 1, where a and /3 label the two 
coexisting bulk phases, good results are obtained. How- 
ever, when we consider a difference in viscosity between 
the two components, = 3 and Tp = 0.7, the approach 
does not work well. In the next section we shall explain 
why not and describe a method to overcome the problem. 

Fig. [1] shows a drop resting on a solid surface. In gen- 
eral, the fluid-solid {jas and ^ps) and fluid-fluid (7) sur- 
face tensions are different. At the contact point A the 
balance of forces is described by Young's law 



cosr 



7 



(20) 



where 9'^'^ defines the equilibrium contact angle. 

We performed simulations to verify this relation. The 
lattice Boltzmann system size was set to 300 x 100 lat- 
tice units and the parameters used were a = 0.04, = 1, 
F = 0.5, At = Ax = 1 , and k ~ 0.04, giving an inter- 
face width of W = 2^2^/0 = 2.8 lattice sites and an 
interfacial tension of 7 = 0.038 in lattice units. 

The relaxation parameter Tp was determined by 



= T/3 + ^ {Ta - Tp) 



(21) 



such that it changed smoothly through the interface and 
had the bulk values Tq and Tp in the two bulk phases. 
The values Tq = 3 and Tp = 0.7 were chosen to give a 
viscosity ratio of Ri, = (tq — 0.5) /{Tp — 0.5) = 12.5. 

Non-slip boundary conditions at the walls were 
achieved using a standard mid-link, bounce-back method 
[2]] | . The contact angle was varied by changing the gradi- 
ent of the order parameter at the solid boundary, 9±(/>|(,. 
Initially, a semi-circle of fluid of radius i? = 35 was placed 
on the surface with a contact angle of 90°. The bound- 
ary conditions were set to i9_l<?I>|j — and the system 
was evolved for 3 x 10^ time-steps, such that it reached 
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equilibrium, d^ff'lb then quasi-statically increased 
over the course of 10^ time steps and the variation in 
the contact angle was measured. The process was re- 
peated, but this time decreasing 9^01^ until the surface 
was completely wet. 

Numerical measurements of contact angle were per- 
formed by matching the interface to the arc of a circle 
[1^. Specifically, each link between neighbouring lattice 
sites was examined to see if at one end it was fluid a 
{(p > 0) and at the other end fluid /?((/>< 0). If this was 
the case then a linear interpolation was used to predict 
the point on the link where </> = 0. These crossing points 
are illustrated by the crosses in Fig. [1] (note that this 
grid is not to scale). The equilibrium contact angle was 
estimated by performing a least squares fit of the cross- 
ing points to a circular section and then calculating the 
contact angle the section made with the surface. 

Results for the equilibrium contact angle for different 
values of Tq, and are shown in Fig. [2] and compared to 
the exact result 



9x01 b = 2sgn r-"- 



cos 



(f){'-(i)}r><-) 



where a = arccos (sin 0'^'?) . The agreement is good for 
= = 1 (Fig. UK a)), but there are large errors for 
Ta 7^ Tj3 (Fig. mb)). In the next section we discuss the 
reasons behind this discrepancy. 



V. THE ORIGINS OF SPURIOUS CURRENTS 
NEAR THE CONTACT POINT 

The reason that the single relaxation time lattice 
Boltzmann approach gives an incorrect equilibrium con- 
tact angle for Tp 7^ 1 is that near to the contact point 
there are strong spurious velocities which continuously 
push the system out of equilibrium and result in the de- 
formation of the drop. Note that even when the surface 
is neutrally wetting, 9j_0|^ = 0, the wetting angle mea- 
sured numerically is ~ 10° larger than the expected value 
of = 90°. We focus on this case and replace the drop 
geometry with a simpler system consisting of a stripe of 
component a between two neutrally wetting walls, as de- 
picted in Fig. [Sja). This simulation was performed for a 
system of size 28 x 15 lattice units using a large relaxation 
parameter Tp = Ta = T/3 = 10. Note that we do observe 
the correct 90° contact angle in this case, but only be- 
cause the viscosities of the two phases are the same. If 
we set Ta > Tp then the stripe becomes barrel shaped, 
consistent with the measurements reported in Fig. HJa). 

The black arrows in Fig. [3Ka) show the steady-state, 
spurious velocity field in the system. The magnitude of 
the velocities is typically of order 10~^c. We have iden- 
tified two contributions to the spurious velocities; one 
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FIG. 2: The equilibrium contact angle as a function of the gra- 
dient in (j> at the boundary (measured in lattice units), d±4>\b-: 
for (a) Ta = Tf) = l.Q and (b) Ta = 3, r/3 = 0.7. Squares 
and crosses are simulation results obtained using single and 
multiple-relaxation-time lattice Boltzmann algorithms respec- 
tively. The solid curve is the theoretical expression p2|l . 



from long range effects, a second from the bounce back 
boundary conditions. We discuss each in turn. 



A. Spurious velocities from long range effects 

For a imbounded, system at steady-state, it is possi- 
ble, by iterating the lattice Boltzmann evolution equa- 
tion (jlOp . to write the particle distribution function at 
any given lattice point in terms of equilibrium distribu- 
tions along lines defined by the lattice velocity vector 
directions: 
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FIG. 3: (a) The spurious flow field in steady state in a system 
consisting of a stripe of one fluid (light region) surrounded by 
another fluid (dark regions) held between two walls. Peri- 
odic boundary conditions are applied in the a--direction. The 
parameters used were Ta = rp = 10. (b) (Color online) An en- 
largement of the region near the interface indicating how the 
particle distribution function at A depends on contributions 
from lattice nodes in the lattice velocity vector directions. 



Note that when Tp = 1 this reduces to /i(r) = f^'^{r — 
AtSi). This special case could be described as 'local', 
in the sense that the particle distribution function only 
depends on the equilibria of its neighbours. 

When Tp > 1, contributions to the sum in Eq. 
decay with a characteristic length 



>i 



In 



Ax. 



(24) 



This diverges as Tp — > oo. When Tp < 1, each term in 



the sum (j23p alternates in sign and the envelope decays 
exponentially with length scale 



<i 



In 



Aa; 



(25) 



which diverges as Tp ^ (This limit makes sense because 
it corresponds to the viscosity in Eq. becoming zero.) 
Therefore, for Tp both high and low, the distribution at 
any given point is dependent on other nodes a long dis- 
tance away. 

Even without the presence of solid boundaries, it has 
been noted that these long range interactions can give rise 



to large spurious currents around curved interfaces [l7|. 
In particular, the size of the spurious velocities scales as 
Tp for large Tp, so choosing Tp large is not advisable. The 
problem becomes even more acute in the presence of solid 
boundaries. 

Figurc[3ljb) schematically shows the situation for a par- 
ticular lattice node A close to the contact point. Because 
we focus on the case of neutrally wetting walls (contact 
angle 90°), this allows us to compare two systems we ex- 
pect to have the same steady state configuration, i.e. a 
rectangular stripe of fluid at rest. In the first, simpler 
system we imagine that the walls have been replaced by 
periodic boundary conditions in the y-direction. In this 
case, contributions to Eq. (|23p radiate out indefinitely in 
all 8 directions. This is illustrated in Fig. ^h) by the 
filled circles, each of which denote a term in the sum on 
the right hand side of Eq. ([23]) . (Note, the arrows indi- 
cate the direction of motion for the distribution function 
from which each term originates.) Because the bound- 
aries have been removed, the symmetry of the system 
implies that any point above or below A at the same 
value of X must behave exactly the same as A. Since 
there can be no net flux of fluid across any given plane 
(as this would lead to non-conservation of mass), then 
in steady state, the system must be at rest. It can be 
numerically confirmed that no steady state spurious ve- 
locities are generated in this case. 

In the second system, the periodic boundary conditions 
are replaced with bounce back boundary conditions, and 
we now discuss what implications this has. In particu- 
lar. Eq. ([23|l needs to be modified to take into account 
the fact that the contributions to the sum do not extend 
indefinitely in all 8 directions because of the presence 
of the boundaries. This is illustrated at point B in Fig. 
|31[b), where the remaining terms in the sum, which would 
have originated in the direction of C, come, in fact, from 
a reflected branch in the direction of A. Because of the 
interface in the system, these two paths give significatly 
different contributions to the sum. Therefore, the reflec- 
tion of branches at the boundary breaks the symmetry of 
the system in the y direction. (This can clearly be seen 
if we compare point A to a point on the boundary, where 
the incoming branches arc reflected immediately.) This 
broken symmetry leads to the generation of the spurious, 
steady state velocities. 

The range around the contact point over which this 
spurious force is active is determined by the decay lengths 
A in Eqs. (|24|) and p5|) . For high or low viscosities A is 
large and this is a long range effect. Correcting for it at 
the boundary involves extrapolating into the surface to 
predict, for example, the density variation down the dot- 
ted branches near to C. This is relatively simple for the 
90° case shown, as the unknown nodes can be obtained by 
reflecting the system across the boundary. However, in 
the case of arbitrary contact angle the situation is much 
more complicated and a general solution is far from clear. 

Thus implementing solid boundaries in the multi- 
component BGK lattice Boltzmann is only advisable in 
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FIG. 4: (Color online) (a) The mid-link bounce-back method, 
before and after the streaming step, along one lattice di- 
rection, (b) The lattice node bounce-back method, (c-f) 
Schematic diagrams showing the particle distribution func- 
tions for (c) a single-component fluid, (d) a single-component 
fluid under shear, (e) particle distribution function in the in- 
terface of a binary fluid before and (f) after bounce back. 



the case To = 1. (For the same reason it is best to choose 
T0 = 1 |22l|.) However, to simulate phases with dilferent 
viscosities it is not possible to use Tp = 1 in both phases. 
In Sec. IVII we discuss how using a multiplc-rclaxation- 
time method can be used to overcome this restriction. 



B. Spurious velocities from bounce-back boundary 
conditions 



streaming the population a moves down the link, is re- 
flected, and then travels back to the node it came from 
but moving in the opposite direction. (Fig. Sfb) shows 
an alternative bounce-back scheme in which the wall lies 
on the boundary nodes. The arguments presented here 
are equally valid for this case.) 

Fig. IHc) shows the distribution for a typical node in 
a system describing a single component fluid at equilib- 
rium. This is a representation of the Maxwell Boltzmann 
distribution using a discrete number of velocity vectors. 
The bounce-back method acts at the boundaries of the 
system to switch some of the populations to travel in the 
opposite direction. Because this distribution is invariant 
under a parity transformation (that is replacing velocity 
vectors v by — v) then the correct distribution is pre- 
served by the bounce-back boundary. When the system 
is in a state of shear the particle distribution function, 
in the rest frame of the fluid, is as depicted in Fig. [D^d). 
While this is clearly not isotropic, it still preserves in- 
variance under the parity transformation and hence the 
bounce-back approach is still valid. (In fact, if non-slip 
boundary conditions are enforced by setting the parti- 
cle distribution at the boundary to its equilibrium value 
at rest, an inaccurate shear velocity profile is obtained. 
This is because the shear-induced distortion of the dis- 
tribution function is not preserved at the boundaries.) 

We now return to the binary fluid case. Fig. UJ^e) shows 
the distribution function for a typical node, lying at rest 
in a fluid-fluid interface. (In this particular example the 
interface lies perpendicular to the surface, hence the dis- 
tribution has an up-down symmetry.) Note that because 
of the K, terms in the pressure tensor ^ the parity invari- 
ance is broken. We consider the case when the position of 
the node is lS.x/2 above a boundary and mid-link bounce- 
back boundary conditions are being employed. Fig. IDJf) 
illustrates the situation after the bounce back step, and it 
clearly shows that the new distribution is not the same as 
in (e) (in particular, vector a is shortened and b length- 
ened). Therefore, bounce-back collisions result in the 
system continually being pushed out of equilibrium, lead- 
ing to the generation of spurious velocities on or near to 
boundaries close to the interface between the two phases. 

For Tp = 1, the distribution function is automatically 
set to its equilibrium value at each time step. While 
the typical distribution function in the interface is not 
invariant under a parity transformation, the equilibrium 
distribution is. Hence the spurious velocities caused by 
bounce back boundary conditions are suppressed in this 
case. 



A second source of spurious velocities at the contact 
point is the bounce-back boundary conditions. To un- 
derstand what goes wrong with the boundary conditions 
for the binary lattice Boltzmann model it is first impor- 
tant to understand why they work well when simulating 
a single-component fluid. Fig. [D^a) illustrates the mid- 
link, bounce-back process for one dimension. During the 



VI. SUPPRESSING THE SPURIOUS 
CURRENTS 

We now describe how an algorithm based on the 
multiple-relaxation-time lattice Boltzmann approach can 
be used to significantly reduce the spurious currents at 
the contact point. The approach comprises the following 
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four steps: 

Step 1: Calculate the density, concentration and ve- 
locity using the moments defined by Eq. 

Step 2: Set the velocity of boundary nodes to zero 
when calculating the equilibrium distribution function. 

Step 3: Use the multiple-relaxation-time lattice Boltz- 
mann method described in Sec. lIIIBl to perform the col- 
lision step. 

Step 4: Perform the streaming step with bounce-back 
at the boundaries. 

To justify this choice we need to consider the hydro- 
dynamic and the non-hydrodynamic modes separately. 
From a hydrodynamic point of view, a non-slip boundary 
fixes the velocity at that boundary to zero. If the fluid is 
incompressible, then near to the boundary the fluid flow 
profile can always be approximated by a shear profile, 
with streamlines parallel to the surface. In section|Vj3 we 
argued that for a shear profile, the hydrodynamic modes, 
which are represented by the difference between Figs, 
m^c) and (d) , are invariant under a parity transformation, 
and so are preserved by the bounce-back method. An- 
other interpretation of this is that, for the hydrodynamic 
modes, the contributions at A in Fig. [3] from the reflected 
branches are the same as those when periodic boundary 
conditions arc considered (provided Ar^ is smaller than 
the length scale over which the shear profile approxima- 
tion breaks down as we move away from the boundary). 
In summary, the hydrodynamic modes behave correctly 
for both the single and multiple-relaxation-time lattice 
Boltzmann models. The difference lies in the treatment 
of the non-hydrodynamic modes. 

In the multiple-relaxation-time scheme, because we 
have chosen Tp = 1 for the non-hydrodynamic modes 
(this choice corresponds to the I's in Eq. p^ ). the spu- 
rious velocities generated from long ranged effects in the 
bulk are immediately removed. This justifies the use of 
the multiple-relaxation-time algorithm in Step 3. The 
only potential problem that remains is on the boundary 
nodes themselves. As discussed at the end of section IVB. 
the distribution function in the interface between fluid 
phases is not invariant under a parity transformation. 
This generates spurious velocities at any boundary node 
in an interface immediately after the streaming step 4. 
This problem is remedied by the introduction of step 2. 
Note that this step is consistent with non-slip boundary 
conditions and docs not affect the hydrodynamic correla- 
tions (i.e. the non symmetric distribution in Fig. IHd)). 

For a system with variable viscosity it would seem nec- 
essary to recalculate the collision matrix C = M~^SM 
in P?|) at each lattice node and at each time-step. This 
would be extremely slow computationally and not very 
practical. The approach we take is to create a lookup 
tabic with ~ 10* different values of viscosity and simply 
pick the closest match. 

We find that implementing multiple-relaxation-time 
lattice Boltzmann with appropriate boundary conditions 
leads to a significant improvement in the accuracy of the 
equilibrium contact angle. The results of simulations for 



Ta = 3, T/3 = 0.7, are denoted by the crosses in Fig. 
dja), and show very good agreement with the theoreti- 
cal, dashed curve. Deviation is only noticeable at small 
contact angles. This happens for two reasons: Firstly, 
the dynamics of drop wetting become very slow as the 
equilibrium contact angle becomes small, and so the as- 
sumption of quasi-static equilibrium, as the gradient in 
(f) is slowly reduced, breaks down. (Tanner's law (23j 
states that for a completely wetting surface in two dimen- 
sions, the size of a drop spreading on the surface scales 
as r ~ i^/*" in time.) Secondly, the finite width of the 
interface W, which is neglected when assuming that the 
drop should be made up of a circular section, becomes 
comparable to the height of the drop. 



VII. A DYNAMICAL TEST: CAPILLARY 
FILLING 

We have shown that the equilibrium contact angle is 
accurately recovered in lattice Boltzmann simulations for 
a binary system with different viscosities only when a 
multiple-relaxation-time approach is employed. While 
this highlights a problem with the single relaxation time 
binary lattice Boltzmann model, contact angle measure- 
ment is a static problem. In this section we concentrate 
on the dynamics of a fluid penetrating a smooth mi- 
croclianncl to measure how the advancing contact angle 
changes with the velocity of the fluid interface. 

Fluid is pulled into a hydrophilic capillary by the 
Laplace pressure across the interface. Balancing this 
against the viscous drag of the fluid column gives Wash- 
burn's law [23| describing the variation of the length / of 
fluid in the capillary with time t 

/^^(^^)(^ + M, (26) 

where h is the width of the capillary, 9" is the advanc- 
ing contact angle and is an integration constant. Note 
that it is appropriate to use, not the static, but the ad- 
vancing contact angle 6*° , as this controls the curvature of 
the interface and hence the Laplace pressure. Moreover, 
Eq. ((26l) assumes that there is no resistance to motion 
from any fluid already in the capillary. 

Numerical results showing capillary filling of a two 
dimensional capillary are presented in Fig. [5] for both 
the single and the multiple-relaxation-time lattice Boltz- 
mann algorithms. The simulations are for a channel of 
length L 640 and width /i = 50. Reservoirs (480 x 200) 
of the two components are attached at each end of the 
capillary. The two reservoirs are connected to ensure 
that they have the same pressure. The other parameters 
of the model are chosen to give 9'^'^ = 60°, 7 = 0.0188, 
and pi ~ p2 = 1-0. Fluid a, with viscosity vi = 0.83 is 
taken to displace fluid /3 with viscosity 1/2 = 0.03. Re- 
sults are shown for mobilities M = 0.05, Af = 0.1 and 
M = 0.5. 
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FIG. 5: Lattice Boltzmann simulation results for capillary filling in smooth channels. Left: the length squared, in lattice 
units, of the column of the filling component (fluid a) plotted against time for (a) single and (c) multiple-relaxation-time 
lattice Boltzmann simulations. The circles are the simulation results and the solid lines are fits to Washburn's law. Right: the 
advancing contact angle of the liquid-liquid interface for (b) single and (d) multiple-relaxation-time simulations. The crosses 
are the simulation results and the solid lines are linear fits of cos 9°" to the capillary number [2^. 

gorithni, the advancing contact angle tends to the cor- 
rect value as Ca — *■ 0. We obtain 9a\ca^o = 58°, 60° 
and 60° for M = 0.05, 0.1 and 0.5 respectively. For 
the BGK method, however, 6''^|ca-.o = 75°, 73° and 
72° for M = 0.05, 0.1 and 0.5, and the advancing an- 
gle is consistently higher for all values of Ca than for 
the multi-relaxation-time solution. (This result agrees 
with that presented in Fig.[2l which shows that the mea- 
sured equilibrium contact angle is too high in the single- 
relaxation-time approach.) Since the speed of capillary 
filling depends on cos^" (HH), capillary filling is consid- 
erably slower using this method. 

In diffuse interface models of binary fluids the contact 
line singularity is relieved by inter-diffusion of the two 
fluid components. This is governed by the mobility M; 
as is apparent from Fig. [5l increasing M increases the 
rate of diffusion and hence the velocity of the contact 
line. Therefore the parameter M can be used to tune the 
effective slip length Is- 



The solid lines in Fig.ISJa) and (c) are flts to the Wash- 
burn law. At later times, once inertia! effects have be- 
come negligible, both the single and multiple-relaxation- 
time lattice Boltzmann simulations satisfy the Washburn 
scaling. However, the two methods provide different 
quantitative results, in particular, capillary flUing is con- 
siderably slower for the single-relaxation-time method. 
The difference can be explained by considering how the 
advancing contact angle varies with the velocity of the 
interface, . 

To lowest order in the capillary number, Ca= v^v/^^ 
the advancing contact angle is related to the equilibrium 
angle and the capillary number by [25| 

cosr = cos^'' - Calog(A'i//,), (27) 

where K \s a. constant, L is the length scale of the sys- 
tem and Is is the effective slip length at the three phase 
contact line. Figs. [Hb) and (d) show the expected lin- 
ear decrease of the measured contact angle with capil- 
lary number [251 . For the multiple-relaxation-time al- 
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VIII. SUMMARY AND CONCLUSIONS 

We have shown that, if the lattice Bohzmann relax- 
ation parameter Tp 7^ 1, strong spurious currents drive 
the contact line out of equilibrium even in a system at 
rest. This means that the algorithm gives incorrect values 
for the contact angle. In many applications it is possible 
to choose Tp = 1, thus avoiding the problem. However 
this parameter controls the fluid viscosity and so cannot 
be held at unity for both phases of a binary fluid if the 
two components have different viscosities. 

We demonstrate that the spurious currents arise pri- 
marily from two effects. The first is a long-range contri- 
bution to the equilibrium distribution function near the 
contact line that effectively originates in the incorrect 
phase. The second is the bounce-back boundary condi- 
tions which drive the interface out of equilibrium. 

Aiming to reduce the unwanted velocities we propose a 
revised lattice Boltzmann method, based on a multiple- 
relaxation-time algorithm. We show that the simulations 
then agree well with the analytical result for the equi- 
librium contact angle. Moreover the dynamic, advanc- 
ing contact angle shows the expected physical behaviour, 
with a slip length that depends on the diffusivity of the 
fluid. 

Using this method it will be possible to perform accu- 
rate simulations of a binary fluid where the two com- 
ponents have different viscosities. Examples of prob- 
lems where the algorithm will prove useful include vis- 
cous fingering, rolling of viscous drops and capillary fill- 
ing. Moreover, it can provide a useful approximation 
to a liquid-gas system in the limit that evaporation- 
condensation is not important. 



based on reducing the magnitude of spurious velocities 
near to interfaces (a detailed account is given in [13 )• 

The equilibrium distributions can be written in the 
form 



3 



pUaUfj 



2r? 



(j)UaUp 1 , 



(Al) 

for I = 1, .., 8, where wi-4 = ^, W5-8 = j^, and summation 
over repeated indices is assumed. Other parameters are 



w 



1-2 

xy 



yy 



3-4 



0, w', 



xy 



Uf_4 



W 



yy _ 



1-2 



5-8 



yy 

W5-8 = -- 



i_4 - u, LC5,6 - 4' and g = - i. 
The i = stationary values are chosen to conserve the 
concentration of each species: 



i=i 



During the lattice Boltzmann procedure, it is necessary 
to calculate numerically both derivatives {e.g. dxP in the 
equilibrium distribution (jAip ) and the Laplacian {e.g. 
to obtain the chemical potential ([3])). These continuous 
quantities are calculated from stencils, discrete operators 
which use neighbouring lattice sites. The best choice of 
stencils to reduce spurious velocities is given by 



APPENDIX A: THE CHOICE EQUILIBRIUM 
DISTRIBUTION 



dx- 



12Ai 



-1 1 
-4 4 
-1 1 



6(A2;) 



1 4 1 
4 -20 4 
1 4 1 



(A3) 



We list the best choice of equilibrium distributions and 
stencils for calculating spatial derivatives for the lattice 
Boltzmann algorithms we have considered in this paper, 



The values in these matrices denote the weights given to 
a particular quantity on a lattice node (the central entry) 
and on the surrounding eight lattice points. 
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